global M;
M = 8;
Fs = 16000;
Wp = 3400/8000; Ws = 4600/8000;
Rp = 0.3; Rs = 60;
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);
                   
T = 1/Fs;                     % Sample time
L = 2048;                     % Length of signal
t = (0:L-1)*T;               
x =  zeros(1,2048); 
x(1) = 2^(M-1);
NFFT = 2^nextpow2(L);

xq = round((2^(M-1)-1)*x) / 2^(M-1);

[num,den] = ellip(n,Rp,Rs,Wp);
y0 = filter(num,den,xq);
FY0 = fft(y0,NFFT)/L;

[h,w] = freqz(num,den,NFFT,Fs);
[z,p,k] = tf2zp(num,den);

zq = round((2^(M-1)-1)*z) / 2^(M-1);
pq = round((2^(M-1)-1)*p) / 2^(M-1);
%pq2 = round((2^(M-1)-1)*p*128) / (2^(M-1)*128)
kq = round((2^(M-1)-1)*k) / 2^(M-1);
%kq2 = round((2^(M-1)-1)*k*128) / (2^(M-1)*128)

[numz,denz] = zp2tf(zq,pq,kq);
numq = round((2^(M-1)-1)*numz) / 2^(M-1);
denq = round((2^(M-1)-1)*denz) / 2^(M-1);
[hq,wq] = freqz(numq,denq,NFFT,Fs);
y = filter(numq,denq,xq);
FY = fft(y,NFFT)/L;


[sos,g] = zp2sos(zq,pq,kq);
[sos1,g1] = zp2sos(zq,pq,kq,'up','inf');  %infinity-norm scaling with up-ordering
[sos2,g2] = zp2sos(zq,pq,kq,'down','two');  %2-norm scaling with down-ordering


gq = round((2^(M-1)-1)*g) / 2^(M-1);
g1q = round((2^(M-1)-1)*g1) / 2^(M-1);
g2q = round((2^(M-1)-1)*g2) / 2^(M-1);

ys = gq*sosfilter(sos(3,1:3), sos(3,4:6),sosfilter(sos(2,1:3),sos(2,4:6),sosfilter(sos(1,1:3),sos(1,4:6),xq)));
FYs = fft(ys, NFFT)/L;

yinf = g1q*sosfilter(sos1(3,1:3), sos1(3,4:6),sosfilter(sos1(2,1:3),sos1(2,4:6),sosfilter(sos1(1,1:3),sos1(1,4:6),xq)));
FYinf = fft(yinf, NFFT)/L;

y2 = g2q*sosfilter(sos2(3,1:3), sos2(3,4:6),sosfilter(sos2(2,1:3),sos2(2,4:6),sosfilter(sos2(1,1:3),sos2(1,4:6),xq)));
FY2 = fft(y2, NFFT)/L;




f = Fs/2*linspace(0,1,NFFT);
hview = [hq,h];

% subplot(2,1,1);
plot(f,20*log(abs(FY)),'--',f,20*log(abs(FYinf)),'--',f,20*log(abs(FY2)),':',f,20*log(abs(FYs)),f,20*log(abs(FY0)));
% grid on;

% subplot(2,1,2);
% semilogx(w*Fs/(2*pi), 20*log10(abs(hview)));
% grid on;
